#=

=#

using LinearAlgebra


"""

"""
function fakedata(Ls)
    hop = zeros(Ls, Ls)
    intup = zeros(Ls, Ls)
    intdn = zeros(Ls, Ls)
    for i in Base.OneTo(Ls-1)
        hop[i, i+1] = -1.0
        hop[i+1, i] = -1.0
    end
    for i in Base.OneTo(Ls-2)
        hop[i, i+2] = -1.0
        hop[i+2, i] = -1.0
    end
    ehk = exp(-0.1*hop)
    #
    U = 1.0
    dt = 0.1
    lam = acosh(exp(U*dt/2))
    aux = 2(round.(rand(Ls)) .- 0.5)
    intup += Diagonal(-aux*lam)
    intdn += Diagonal(aux*lam)
    #println(intup)
    #println(intdn)
    intup = exp(intup)
    intdn = exp(intdn)
    return intup*ehk, intdn*ehk
end


"""
size(B) = [2*Nx,2*Nx]
"""
function bmat_Ising(hk, cfg, dt; U=1.0)
    #ehk = [0 1; 1 0]
    ehk = exp(-dt*hk)
    #
    ehk2 = kron([1 0; 0 1], ehk)
    #
    #U = 1.0
    lam = acosh(exp(U*dt/2))
    intup = -cfg*lam
    intdn = cfg*lam
    eV = Diagonal(vcat(intup, intdn))
    eV = exp(eV)
    return eV * ehk2
end


"""
实部虚部分开处理，
∂L/∂nx = (∂L/∂B ∂B/∂Br) ∂Br/∂nx + (∂L/∂B ∂B/∂Bi) ∂Bi/∂nx = ∂L/∂B ∂Br/∂nx + i ∂L/∂B ∂Bi/∂nx
上式的实部是 real(∂L/∂B) * (∂Br/∂nx)  - imag(∂L/∂B) * (∂Bi/∂nx)
"""
function bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)
    #
    #ehk = [0 1.0+1.0im; 1.0-1.0im 0]
    #ehk = exp(-dt*hk)
    #if Z2
    #    ehk2 = kron([1 0; 0 1], ehk)
    #else
    #    ehk2 = ehk
    #end
    #这里的错误需要修改
    #exp([nz*lam nx*lam-ny*lam; nx*lam+ny*lam -nz*lam])
    #不应该是Ui，应该是每个格点的lam
    lams = acosh.(exp.(Ui*dt/2))
    clam = -cfg.*lams
    #
    M11 = @. exp(clam)*(1+nz)/2 - exp(-clam)*(nz-1)/2
    M12 = @. (nx-ny*im)*exp(clam)/2 - (nx-ny*im)*exp(-clam)/2
    M21 = @. (nx+ny*im)*exp(clam)/2 - (nx+ny*im)*exp(-clam)/2
    M22 = @. exp(clam)*(1-nz)/2 - exp(-clam)*(-1-nz)/2
    #println(typeof(cfg))
    #println(exp.(Ui*dt/2)) 
    #ints = kron([1 0; 0 0], Diagonal(exp.(Uz)))
    #coef = Diagonal(cfg.*lams)
    ints = kron([1 0; 0 0], Diagonal(M11))
    ints += kron([0 0; 0 1], Diagonal(M22))
    ints += kron([0 1; 0 0], Diagonal(M12))
    ints += kron([0 0; 1 0], Diagonal(M21))
    bmat = ints * ehk
    return real(bmat), imag(bmat)
end


"""
每一个位置分开算，加速计算的过程
"""
function bmat_IsingADX(ehk, ::Val{CFG}, Ui, nx, ny, nz, dt) where CFG
    lams = acosh(exp(Ui*dt/2))
    clam = -CFG*lams
    #
    M11 = exp(clam)*(1+nz)/2 - exp(-clam)*(nz-1)/2
    M12 = (nx-ny*im)*exp(clam)/2 - (nx-ny*im)*exp(-clam)/2
    M21 = (nx+ny*im)*exp(clam)/2 - (nx+ny*im)*exp(-clam)/2
    M22 = exp(clam)*(1-nz)/2 - exp(-clam)*(-1-nz)/2
    #
    ret = [M11 M12; M21 M22] * ehk
    return vcat(real(ret), imag(ret))#
end


"""
计算每个位置但是所有虚时
"""
function bmat_IsingADTX(ehk, ths, Ui, nx, ny, nz, dt)
    Nt = length(ths)
    DNx = size(ehk)[2]
    lams = acosh(exp(Ui*dt/2))
    clam = -ths*lams
    #
    M11 = @. exp(clam)*(1+nz)/2 - exp(-clam)*(nz-1)/2
    M12 = @. (nx-ny*im)*exp(clam)/2 - (nx-ny*im)*exp(-clam)/2
    raw1 = hcat(M11, M12)
    raw1 = raw1 * ehk
    #raw1 = reshape(raw1, 1, Nt, DNx)
    M21 = @. (nx+ny*im)*exp(clam)/2 - (nx+ny*im)*exp(-clam)/2
    M22 = @. exp(clam)*(1-nz)/2 - exp(-clam)*(-1-nz)/2
    raw2 = hcat(M21, M22)
    raw2 = raw2 * ehk
    #raw2 = reshape(raw2, 1, Nt, DNx)
    #
    return vcat(real(raw1), real(raw2), imag(raw1), imag(raw2))#cat(real(raw1), real(raw2), imag(raw1), imag(raw2); dims=2)
end


"""
需要ehk
"""
function bmat_IsingND(ehk, cfg, Ui, nx, ny, nz, dt)
    bmr, bmi = bmat_IsingAD(ehk, cfg, Ui, nx, ny, nz, dt)
    return complex.(bmr, bmi)
end


"""
Δ = (eV'e-V - 1)
"""
function Δmat_Ising(Nx, ni, ncfg, ocfg, dt; U=1.0)
    lam = acosh(exp(U*dt/2))
    del = Diagonal(zeros(Nx*2))
    del[ni, ni] = exp(-(ncfg-ocfg)*lam) - 1.0
    del[ni+Nx, ni+Nx] = exp((ncfg-ocfg)*lam) - 1.0
    return del
end


function Δmat_IsingND(Nx, ni, ncfg, ocfg, U, nx, ny, nz, dt)
    #
    lams = acosh.(exp.(U*dt/2))
    clam = -(ncfg-ocfg)*lams
    delc = zeros(ComplexF64, 2, 2)
    #
    #M11 = clam*nz
    #M12 = clam*nx - clam*ny*im
    #M21 = adjoint(M12)
    #M22 = -clam*nz
    #delcompact = exp([M11 M12; M21 M22]) - I(2)
    delc[1, 1] = exp(clam)*(1+nz)/2 - exp(-clam)*(nz-1)/2 - 1.0
    delc[1, 2] = (nx-ny*im)*exp(clam)/2 - (nx-ny*im)*exp(-clam)/2
    delc[2, 1] = (nx+ny*im)*exp(clam)/2 - (nx+ny*im)*exp(-clam)/2
    delc[2, 2] = exp(clam)*(1-nz)/2 - exp(-clam)*(-1-nz)/2 - 1.0
    #delcompact = [M11 M12; M21 M22] - I(2)
    #
    #del[ni, ni] = delcompact[1, 1]
    #del[ni, ni+Nx] = delcompact[1, 2]
    #del[ni+Nx, ni] = delcompact[2, 1]
    #del[ni+Nx, ni+Nx] = delcompact[2, 2]
    return [ni, ni+Nx], delc
end


"""
有多余规范的HS
``````
e^{-τU(n↑n↓-0.5n↑-0.5↓)} = 0.5 y1^{an↑+bn↓} + 0.5 y2^{-cn↑-dn↓}
y1 = e^{τU/2} y2 = e^{-τU/2}
y1^{a} =  y1a
y2^{-c} =  y1^{c} = 2y1 - y1a
y1^{b} = (1 - 2*y1^2 + y1a*y1)/(y1a - y1)
y2^{-d} = y1^{d} = (y1a*y1 - 1)/(y1a - y1)
#指向-1的时候，取a和b，+1取c和d
#b控制的概率还需要在其他位置进行更新
"""
function bmat_Gauge1ADX(ehk, ::Val{-1}, Ui, ra, ia, rb, ib, dt)
    #lams = acosh(exp(Ui*dt/2))
    y1 = exp(dt*Ui/2)
    y1a = y1^(ra + ia*im)
    pval = rb + ib*im
    #
    M11 = y1a
    M12 = 0.0
    M21 = 0.0
    M22 = (1 - pval - y1^2 + pval*y1a*y1)/(pval*y1a - pval*y1)
    #
    ret = [M11 M12; M21 M22] * ehk
    return vcat(real(ret), imag(ret))
end


"""
有多余规范的HS
``````
e^{-τU(n↑n↓-0.5n↑-0.5↓)} = 0.5 y1^{an↑+bn↓} + 0.5 y2^{-cn↑-dn↓}
y1 = e^{τU/2} y2 = e^{-τU/2}
y1^{a} =  y1a
y2^{-c} =  y1^{c} = 2y1 - y1a
y1^{b} = (1 - 2*y1^2 + y1a*y1)/(y1a - y1)
y2^{-d} = y1^{d} = (y1a*y1 - 1)/(y1a - y1)
#指向-1的时候，取a和b，+1取c和d
"""
function bmat_Gauge1ADX(ehk, ::Val{1}, Ui, ra, ia, rb, ib, dt)
    #lams = acosh(exp(Ui*dt/2))
    y1 = exp(dt*Ui/2)
    y1a = y1^(ra + ia*im)
    pval = rb + ib*im
    #
    M11 = (y1 - pval*y1a) / (1-pval)
    M12 = 0.0
    M21 = 0.0
    M22 = (y1a*y1 - 1)/(y1a - y1)
    #
    ret = [M11 M12; M21 M22] * ehk
    return vcat(real(ret), imag(ret))#
end


"""
有多余规范的HS
``````
e^{-τU(n↑n↓-0.5n↑-0.5↓)} = 0.5 y1^{an↑+bn↓} + 0.5 y2^{-cn↑-dn↓}
y1 = e^{τU/2} y2 = e^{-τU/2}
y1^{a} =  y1a
y2^{-c} =  y1^{c} = 2y1 - y1a
y1^{b} = (1 - 2*y1^2 + y1a*y1)/(y1a - y1)
y2^{-d} = y1^{d} = (y1a*y1 - 1)/(y1a - y1)
#指向-1的时候，取a和b，+1取c和d
"""
function bmat_Gauge1AD(ehk, cfg, Ui, ra, ia, rb, ib, dt)
    #lams = acosh(exp(Ui*dt/2))
    y1 = exp.(dt*Ui/2)
    y1a = @. y1^(ra + ia*im)
    pval = @. rb + ib*im
    #当hs=-1，取y1a, hs=+1取(y1 - pval*y1a) / (1-pval)
    #可以表示成0.5*(hs+1)(2y1 - y1a)-0.5*(hs-1)*y1a
    M11 = @. 0.5*(cfg+1)*(y1 - pval*y1a) / (1-pval) - 0.5*(cfg-1)*y1a
    #当hs=-1，取(1 - 2*y1^2 + y1a*y1)/(y1a - y1)
    #当hs=+1，取(y1a*y1 - 1)/(y1a - y1)
    #可以表示成0.5*(hs+1)(y1a*y1 - 1)/(y1a - y1)-0.5*(hs-1)*(1 - 2*y1^2 + y1a*y1)/(y1a - y1)
    M22 = @. 0.5*(cfg+1)*(y1a*y1 - 1)/(y1a - y1) -
    0.5*(cfg-1)*(1 - pval - y1^2 + pval*y1a*y1)/(pval*y1a - pval*y1)
    #
    ints = kron([1 0; 0 0], Diagonal(M11))
    ints += kron([0 0; 0 1], Diagonal(M22))
    bmat = ints * ehk
    return real(bmat), imag(bmat)
end


"""
有多余规范的HS
"""
function bmat_Gauge1ND(ehk, cfg, Ui, ra, ia, rb, ib, dt)
    bmr, bmi = bmat_Gauge1AD(ehk, cfg, Ui, ra, ia, rb, ib, dt)
    return complex.(bmr, bmi)
end


"""
Δ = (eV'e-V - 1)
"""
function Δmat_Gauge1ND(Nx, ni, ncfg, ocfg, Ui, ra, ia, rb, ib, dt)
    #
    y1 = exp(dt*Ui/2)
    y1a::ComplexF64 = y1^(ra + ia*im)
    pval = rb + ib*im
    #
    delc = zeros(ComplexF64, 2, 2)
    #当hs=-1，取y1a, hs=+1取2y1 - y1a
    #从ocfg=-1，ncfg=1时
    #delc[1,1]=(y1 - pval*y1a)/(1-pval)/y1a - 1.0
    #delc[2,2]=(y1a*y1 - 1)/(y1a - y1)/(1 - pval - y1^2 + pval*y1a*y1)*(pval*y1a - pval*y1)-1.0
    #从ocfg=1，ncfg=-1时
    #delc[1,1]=y1a/(y1 - pval*y1a)*(1-pval) - 1.0
    #delc[2,2]=(1 - pval - y1^2 + pval*y1a*y1)/(pval*y1a - pval*y1)/(y1a*y1 - 1)*(y1a - y1)-1.0
    if (ncfg - ocfg) == 2
        delc[1, 1] = (y1 - pval*y1a)/(1-pval)/y1a - 1.0
        delc[2, 2] = (y1a*y1 - 1)/(y1a - y1)/(1 - pval - y1^2 + pval*y1a*y1)*(pval*y1a - pval*y1)-1.0
        delw = (1-pval) / pval
    else
        delc[1, 1] = y1a/(y1 - pval*y1a)*(1-pval) - 1.0
        delc[2, 2] = (1 - pval - y1^2 + pval*y1a*y1)/(pval*y1a - pval*y1)/(y1a*y1 - 1)*(y1a - y1)-1.0
        delw = pval / (1-pval)
    end
    #
    return [ni, ni+Nx], delc, delw
end


"""
4分量的hs变换
``````
e^{-aA^2} = 1/4 ∑_{l=±1,±2} γ(l) e^{sqrt(-a) η(l) A}
γ(±1) = 1+(√6)/3  η(±1) = ±√(6-2√6)
γ(±2) = 1-(√6)/3  η(±2) = ±√(6+2√6)
这里只处理这一种情况
a = -Ui/2
A = (n↑ - n↓ + ϕ)^2
"""
function _bmat_Quad1ADX(ehk, γ, η, Ui, phi, dt)
    #γ和phi暂时用不到，最后的权重是γ*e^{sqrt(-a) η phi}，这部分手动求导
    sqrtma = sqrt(dt*Ui/2)
    M11 = exp(sqrtma*η)
    M12 = 0.0
    M21 = 0.0
    M22 = exp(-sqrtma*η)
    #
    ret = [M11 M12; M21 M22] * ehk
    return vcat(real(ret), imag(ret))#
end


"""
4分量的hs变换
这里只处理这一种情况
a = -Ui/2
A = (n↑ - n↓ + ϕ)^2
"""
@generated function bmat_Quad1ADX(ehk, ::Val{S}, Ui, phi, dt) where S
    if S === -2
        #1-sqrt(6)/3, -√(6+2√6)
        :(_bmat_Quad1ADX(ehk, 0.18350341907227408, -3.301360247771569, Ui, phi, dt))
    elseif S === -1
        #1+sqrt(6)/3, -√(6-2√6)
        :(_bmat_Quad1ADX(ehk, 1.8164965809277258, -1.049295246550581, Ui, phi, dt))
    elseif S === 1
        #1+sqrt(6)/3, √(6-2√6)
        :(_bmat_Quad1ADX(ehk, 1.8164965809277258, 1.049295246550581, Ui, phi, dt))
    elseif S === 2
        #1-sqrt(6)/3, √(6+2√6)
        :(_bmat_Quad1ADX(ehk, 0.18350341907227408, 3.301360247771569, Ui, phi, dt))
    else
        :(error("value error"))
    end
end


"""
4分量的hs变换
这里只处理这一种情况
a = -Ui/2
A = (n↑ - n↓ + ϕ)^2
最后需要补偿一个化学势项
Uϕn↑ - Uϕn↓
"""
function bmat_Quad1AD(ehk, cfg, Ui, phi, dt)
    #lams = acosh(exp(Ui*dt/2))
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    sqrtma = sqrt.(dt*Ui/2)
    #当hs=-1，取y1a, hs=+1取2y1 - y1a
    #可以表示成0.5*(hs+1)(2y1 - y1a)-0.5*(hs-1)*y1a
    M11 = @. exp(sqrtma*map(x -> ηdict[x], cfg))
    #当hs=-1，取(1 - 2*y1^2 + y1a*y1)/(y1a - y1)
    #当hs=+1，取(y1a*y1 - 1)/(y1a - y1)
    #可以表示成0.5*(hs+1)(y1a*y1 - 1)/(y1a - y1)-0.5*(hs-1)*(1 - 2*y1^2 + y1a*y1)/(y1a - y1)
    M22 = @. exp(-sqrtma*map(x -> ηdict[x], cfg))
    #
    ints = kron([1 0; 0 0], Diagonal(M11))
    ints += kron([0 0; 0 1], Diagonal(M22))
    bmat = ints * ehk
    return real(bmat), imag(bmat)
end


"""
4分量的hs变换
这里只处理这一种情况
a = -Ui/2
A = (n↑ - n↓ + ϕ)^2
"""
function bmat_Quad1ND(ehk, cfg, Ui, phi, dt)
    bmr, bmi = bmat_Quad1AD(ehk, cfg, Ui, phi, dt)
    return complex.(bmr, bmi)
end


"""
Δ = (eV'e-V - 1)
"""
function Δmat_Quad1ND(Nx, ni, ncfg, ocfg, Ui, phi, dt)
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    #
    ηdiff = ηdict[ncfg] - ηdict[ocfg]
    sqrtma = sqrt(dt*Ui/2)
    #
    delc = zeros(ComplexF64, 2, 2)
    delc[1, 1] = exp(sqrtma*ηdiff) - 1.0
    delc[2, 2] = exp(-sqrtma*ηdiff) - 1.0
    #
    γdict = Dict(
        -2 => 0.18350341907227408, 
        -1 => 1.8164965809277258,
        1 => 1.8164965809277258,
        2 => 0.18350341907227408
    )
    delw = exp(sqrtma*ηdiff*phi) * γdict[ncfg] / γdict[ocfg]
    return [ni, ni+Nx], delc, delw
end


"""
4分量的hs变换
``````
e^{-aA^2} = 1/4 ∑_{l=±1,±2} γ(l) e^{sqrt(-a) η(l) A}
γ(±1) = 1+(√6)/3  η(±1) = ±√(6-2√6)
γ(±2) = 1-(√6)/3  η(±2) = ±√(6+2√6)
这里只处理这一种情况
a = Ui/2
A = (n↑ + n↓ + ϕ)^2
"""
function _bmat_Quad2ADX(ehk, γ, η, Ui, phi, dt)
    #γ和phi暂时用不到，最后的权重是γ*e^{sqrt(-a) η phi}，这部分手动求导
    sqrtma = sqrt(complex(-dt*Ui/2, 0.0))
    M11 = exp(sqrtma*η)
    M12 = 0.0
    M21 = 0.0
    M22 = exp(sqrtma*η)
    #
    ret = [M11 M12; M21 M22] * ehk
    return vcat(real(ret), imag(ret))#
end


"""
4分量的hs变换
这里只处理这一种情况
a = Ui/2
A = (n↑ + n↓ -1 + ϕ)^2
"""
@generated function bmat_Quad2ADX(ehk, ::Val{S}, Ui, phi, dt) where S
    if S === -2
        #1-sqrt(6)/3, -√(6+2√6)
        :(_bmat_Quad2ADX(ehk, 0.18350341907227408, -3.301360247771569, Ui, phi, dt))
    elseif S === -1
        #1+sqrt(6)/3, -√(6-2√6)
        :(_bmat_Quad2ADX(ehk, 1.8164965809277258, -1.049295246550581, Ui, phi, dt))
    elseif S === 1
        #1+sqrt(6)/3, √(6-2√6)
        :(_bmat_Quad2ADX(ehk, 1.8164965809277258, 1.049295246550581, Ui, phi, dt))
    elseif S === 2
        #1-sqrt(6)/3, √(6+2√6)
        :(_bmat_Quad2ADX(ehk, 0.18350341907227408, 3.301360247771569, Ui, phi, dt))
    else
        :(error("value error"))
    end
end


"""
4分量的hs变换
这里只处理这一种情况
a = Ui/2
A = (n↑ + n↓ -1 + ϕ)^2
最后需要补偿一个化学势项
-Uϕn↑ -Uϕn↓
"""
function bmat_Quad2AD(ehk, cfg, Ui, phi, dt)
    #lams = acosh(exp(Ui*dt/2))
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    sqrtma = sqrt.(complex.(-dt*Ui/2, 0.0))
    #当hs=-1，取y1a, hs=+1取2y1 - y1a
    #可以表示成0.5*(hs+1)(2y1 - y1a)-0.5*(hs-1)*y1a
    M11 = @. exp(sqrtma*map(x -> ηdict[x], cfg))
    #当hs=-1，取(1 - 2*y1^2 + y1a*y1)/(y1a - y1)
    #当hs=+1，取(y1a*y1 - 1)/(y1a - y1)
    #可以表示成0.5*(hs+1)(y1a*y1 - 1)/(y1a - y1)-0.5*(hs-1)*(1 - 2*y1^2 + y1a*y1)/(y1a - y1)
    M22 = @. exp(sqrtma*map(x -> ηdict[x], cfg))
    #
    ints = kron([1 0; 0 0], Diagonal(M11))
    ints += kron([0 0; 0 1], Diagonal(M22))
    bmat = ints * ehk
    return real(bmat), imag(bmat)
end


"""
4分量的hs变换
这里只处理这一种情况
a = Ui/2
A = (n↑ + n↓ -1 + ϕ)^2
"""
function bmat_Quad2ND(ehk, cfg, Ui, phi, dt)
    bmr, bmi = bmat_Quad2AD(ehk, cfg, Ui, phi, dt)
    return complex.(bmr, bmi)
end


"""
Δ = (eV'e-V - 1)
"""
function Δmat_Quad2ND(Nx, ni, ncfg, ocfg, Ui, phi, dt)
    ηdict = Dict(
        -2 => -3.301360247771569, 
        -1 => -1.049295246550581,
        1 => 1.049295246550581,
        2 => 3.301360247771569
    )
    #
    ηdiff = ηdict[ncfg] - ηdict[ocfg]
    sqrtma = sqrt(complex(-dt*Ui/2, 0.0))
    #
    delc = zeros(ComplexF64, 2, 2)
    delc[1, 1] = exp(sqrtma*ηdiff) - 1.0
    delc[2, 2] = exp(sqrtma*ηdiff) - 1.0
    #
    γdict = Dict(
        -2 => 0.18350341907227408, 
        -1 => 1.8164965809277258,
        1 => 1.8164965809277258,
        2 => 0.18350341907227408
    )
    delw = exp(sqrtma*ηdiff*(phi-1.0)) * γdict[ncfg] / γdict[ocfg]
    return [ni, ni+Nx], delc, delw
end
